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Abstract 

We consider a Levy flyer of order a that starts from a point xq on an interval 
[O, L] with absorbing boundaries. We find a closed-form expression for the 
average number of flights the flyer takes and the total length of the flights 
it travels before it is absorbed. These two quantities are equivalent to the 
mean first passage times for Levy flights and Levy walks, respectively. Using 
fractional differential equations with a Riesz kernel, we find exact analytical 
expressions for both quantities in the continuous limit. We show that nu- 
merical solutions for the discrete Levy processes converge to the continuous 
approximations in all cases except the case of a — > 2, and the cases of xq —* 
and xq —* L. For a > 2, when the second moment of the flight length distri- 
bution exists, our result is replaced by known results of classical diffusion. We 
show that if xq is placed in the vicinity of absorbing boundaries, the average 
total length has a minimum at a = 1, corresponding to the Cauchy distribu- 
tion. We discuss the relevance of this result to the problem of foraging, which 
has received recent attention in the statistical physics literature. 

I. INTRODUCTION 

In the past two decades, Levy flights and Levy walks (see Pl-fTPI) found numerous appli- 
cations in natural sciences. Realizations of Levy nights in physical phenomena are very 
diverse, including fluid dynamics, dynamical systems, and statistical mechanics. 

In general, Levy flights and Levy walks model anomalous diffusion, which is governed 
by rare but extremely large jumps of diffusing particles. Both Levy walks and Levy flights 
are characterized by broad distributions of their step lengths, for which the second moment 
does not exist. Levy walks and Levy flights of order a < 2 have distributions of step lengths 



1 



with diverging moments of order m > a and converging moments of order m < a. Hence, 
the classical central limit theorem, which governs the behavior of the Brownian motion, is 
not applicable. According to the generalized central limit theorem [H,ll|, the probability 



density ^(x, n) of the displacement x of Levy flights converges after many steps to the Levy 
stable distribution of order a: 

1 r°° 

ty(x, n) = — / exp(— nKq a ) cos(qx)dq, (1) 
n Jo 

where £q is the characteristic width of the distribution of a single step and n is number of 
steps. This distribution is a generalization of the Gaussian distribution, and is character- 
ized for asymptotically large displacements by the power law decay of its density with the 
exponent 

fi = a + l. (2) 

There are several definitions of Levy walks and Levy flights, which differ in terms of their 
spatial-temporal correlations (see e.g. flTOf). Here we will restrict ourselves to the definition 
0. Accordingly, we assume that for Levy flights the duration of each step is constant, so 
that velocity is proportional to the step length. Hence, the time of travel is proportional 
to the number of steps. Consequently, for Levy flights, the mean-square displacement does 
not exist as a function of time. This property impedes direct applications of Levy flights to 
physical phenomena. 

In Levy walks, walkers travel with constant velocity, which is independent of the step 
length. Hence, the time of travel is proportional to the total path length. Consequently, 
the mean-square displacement exists as function of time, but grows faster than linearly. 
This property makes Levy walks applicable for modeling superdiffusion. However, the time 
evolution of Levy flights is simpler than that of Levy walks. Hence, in the following we 
will derive our results for Levy flights, keeping in mind that the total path length of the 
Levy flights corresponds to total time of travel in Levy walks. In the continuous limit, 
Levy flight process is described by the superdiffusion equation, which includes integer first 
order derivative with respect to time and fractional Riesz operator with respect to spatial 
coordinates H,fj,|l2T]l6f . Here we will restrict our study to only this class of equations. 



Note that usually anomalous diffusion is modeled by the equations of the Schneider- 
Wyss type HT7[| , which include fractional derivatives with respect to time and usual Laplace 



or Fokker-Planck operators with respect to spatial coordinates P,f^,|TE|-f2~3|] . In this case, the 
presence of absorbing boundaries can be treated the same way as in normal diffusion, since 
after the separation of variables, the solution can be expressed as series of the usual eigen- 



functions of the boundary problem for the Laplace or Fokker-Planck operator [16,21-23 . 

In the absence of boundaries, the generalized central limit theorem allows us to treat 
Levy flight diffusion with the help of fractional differential equations P,^,|6]j9|JT2]--[15| . In the 



presence of boundaries, the validity of fractional derivative formalism is less clear. Note that 
the problem of the discrete Levy flights is finite, since it involves the characteristic width £q 
of the distribution of discrete steps. Thus the problem of Levy flights in the finite domain 
of linear size L must depend on the ratio M = L/£ . The transition from discrete Levy 
flights to the fractional differential equation involves transition £ —>■ 0. Consequently, the 
total number of flights and total path length diverge as powers of M. Since the same is true 
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for the total path length of the Brownian walker, this problem can be solved by introducing 
fractional diffusion coefficient, the same way it is done in the usual diffusion equation. We 
will address this point in Section IV. 

Levy flights in a slab geometry with absorbing boundaries have been used to model the 
transmission of light through cloudy atmosphere ||26|| . Using heuristic arguments confirmed 
by numerical simulations, ref. p6| found the scaling behavior of the transmission probability 



of a photon through a slab of width L and the total geometrical path length of transmitted 
and reflected light. This behavior was experimentally observed in We analytically 

derive an exact expression for the transmission probability in Section IV. 

Very recently ||16|| , the approximate expressions for the mean first passage time for both 
Schneider- Wyss sub-diffusion equation |L7] and superdiffusion equation |H| have been ob- 
tained by separation of variables. The latter case exactly corresponds to the Levy flight 
problem, which we treat here (see Section IV). The problem of Levy walks on a finite inter- 
val with absorbing boundaries has already been addressed in . In that paper, the authors 
used an integral equation approach and performed the Laplace transform in the temporal 
domain. They found the asymptotic behavior of the survival probability, which is related to 
the asymptotic behavior of the first passage time. An alternative approach to the Levy walk 
problem which employs the fractional Kramers equations can be found in JI5 |. However, as 
far as we know, exact expressions for the mean first passage time for Levy flights and Levy 
walks have not been derived. 

The structure of the paper is the following. In Sections II and III, we find the mean 
time of travel before absorption for both discrete Levy flights and Levy walks as solutions 
of Fredholm integral equations of the second kind p8| with a power-law kernel truncated by 
a cut-off at small distance £q. 

In Section IV, we treat these equations in the limit £q — > and reduce our problem to the 
solution of fractional differential equations with Riesz kernels p|, p!2|j2^ , [29[] . These equations 
were previously applied in the plane contact problem |51J of linear creep theory and were 
solved using spectral relationships with Jacoby polynomials p4| , |3Tf and by the Sonin inver- 
sion formula |H],|3^]. We also show how our method is related to the method of separation 
of variables first applied to the partial fractional differential equation of superdiffusion in 
0. 

In Section V, we compare our analytical solutions obtained in the continuous limit C.q — > 
with the numerical solutions of the Fredholm equations obtained for discrete Levy flights. 
We show that fractional differential equations can serve as good approximations for Levy 
flights with absorbing boundaries for a < 2. We also show that these approximations break 
down when a — > 2 and in the vicinity of the absorbing boundaries. 

Finally, in Section VI, we discuss the relevance of our results to the problem of biological 
foraging. Recently, Levy flights have been used to model animal foraging |7|,|3"3"H3"7f and cell 
diffusion P5| . It has been suggested |36| that Levy flights with a = 1 provide the optimal 
strategy of foraging in case of sparse food sites, if any food site can be revisited. This 
suggestion was based on the optimization of foraging efficiency, defined to be the inverse 
of the average total path length of the flyer before the flyer is absorbed by traps randomly 
distributed with certain density in <i-dimensional space. The average total path length has 
been approximated [0 as a product of the average length of a single flight and the average 
number of flights before the flyer is absorbed by traps. It has been shown PB[ that this 
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product has a maximum at a = 1, if the starting point of the flyer xq is selected in the 
vicinity of the absorbing boundary. Here we confirm this result in the one- dimensional case 
using both an analytical expression for the average total length of flights obtained in the 
continuous limit and the numerical solution for the discrete Levy process. We show that for 
the case of xq in the vicinity of the absorbing boundary, discrete and continuous solutions 
have the same power law asymptotic behavior, but their amplitudes are different. As a 
consequence, the continuous limit approximation has an additional minimum at a — > 2, 
which is absent in the discrete case. This finding indicates that the fractional differential 
equation approach to Levy flights breaks down in the vicinity of the absorbing boundary. 

In Appendix A, we derive the fractional differential operator for the Levy flight problem 
with absorbing boundaries. In Appendix B, we derive exact analytical expressions for the 
number of flights and the total length traveled before absorption, using the Sonin inversion 
formula for the Riesz fractional equation. 



II. MEAN NUMBER OF FLIGHTS 

Consider a Levy flight that starts at point xo of the interval [0, L] with absorbing bound- 
aries. The flyer makes independent subsequent flights of variable random lengths i with 
equal probability in both directions. The length of each flight is taken from the power law 
distribution 

P(\£\>r) = (£ /rr, (3) 

where the exponent a can vary between and 2, and is the minimal flight length, which 
serves as lower cutoff of the distribution. The probability density of the flight length is given 
by 

pW ^« (4) 

where is a normalization constant and 8(x) = 1 for x > or otherwise. The 



exponent /i of refs. p5| , |36| is identical to a + 1. When a > 2, the second moment of the 
flight distribution converges and the process becomes equivalent to normal diffusion. As 
soon as the flyer lands outside the interval [0, L], the process is terminated. Instead of the 
probability density Eq. (f|) one could use any power-law decaying density [|,^| regularized at 
small distances £o, including Levy stable distribution Eq. ([[]) with n = 7r/[2r(a) sin(7ra/2)], 
where T(a) = Jo°° t a ~ l e~ t dt is Euler T-function []39[| . This value of n is selected so that the 
asymptotic behavior of Eq. (ED) coincides with Eq. (H) M. We use a truncated power-law 
density for the simplicity of analytical treatment. 

We are interested in two quantities: the average number of flights before absorption (n) 
and the average total distance (S) traveled before absorption. Note that we consider the 
length of the last flight to be equal to the distance from the previous landing point to the 
boundary of the interval which flyer crosses during its last flight. This condition makes our 



problem equivalent to the problem of Levy walks |25j with the time defined to be equal to 
the sum of the flight lengths. 
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Suppose that the probability density of finding a Levy flyer at point x after n flights is 
P n (x). Then the probability density after n + 1 flights is given by the convolution of the 
probability density P n (x) and the probability density of the next flight p(x) given by Eq. (f|) 



P n +i(y) 



p(y - x)P n (x)dx. 



(5) 



Let C a {£o) be an integral operator with kernel p(x — y) which is defined on a function 
f(x) of an interval [O, L] as 



[C a f](y) = 



ot^o f L f{x)6{\y-x\-£ )dx 



\y-x 



a+l 



(6) 



One can see, that C a is a self-conjugate operator with respect to a scalar product (f,g) = 
J L f(x)g(x)dx. It can be shown that for any continuous function /, J L \£ a f\dx < [1 — 
(2/M)°] Jq \f\dx, where M = L/£ . The value [1 - (2/M) a ] can be regarded as the norm of 
the operator C a . It has a physical meaning of the survival probability, i.e. the probability 
for the flight that starts at the center of the interval to stay unabsorbed, which is less than 
one. 

The distribution after n flights is given by 



P n (x) = [C n a P }(x). 



(7) 



The initial probability density of the flyer located at position xq is the Dirac delta function, 
Po{x) = S(x — Xq). The probability that the flyer remains unabsorbed after n flights, 



[C n a P ]{x)dx < [l-(2/M) a Y 



(8) 



decays exponentially with n. The probability that the flyer is absorbed exactly on the n th 
flight is 



Pr, 



o 



and thus, the average number of flights spent by the flyer on the interval is 



(n) = E p n 



n 



71=1 



L 00 
n=0 







(£ a - J) _1 Pol {x)dx. 



(9) 



(10) 



The infinite sum in Eq. (|l0j) converges, since the norm of C a is less than one. Here {C a — 
is the inverse operator with respect to the operator C a — I, and I is the unity operator. 
Since the operator (C a — l)^ 1 is also a self-conjugate operator, and Pq{x) = 8(x — Xq), we 
have 



(n) = (C a -1)- L h (x ) 



(11) 



where h(x) = —1 is the constant function. This equation can be rewritten as a Fredholm 
integral equation of the second kind |28] with the kernel p(xq — X\): 
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(n(x ))=l+[ (n(xi))p(x -xi)dxi. (12) 
Jo 

Equation ( |T2"D can be interpreted as a recursive method to determine (n(xo)). Indeed, the 
average total number of flights for the process that start at x is equal to one (contribution 
from the very first flight, which always takes place) plus the convolution of the average total 
number of flights for the processes that start from all possible landing points X\ of the first 
flight inside the interval and the probability density p(x — X\) to land at these points after 
the first flight. 

In general, consider a quantity Q(x ) = (Y^LiQi)) where % = g(xj_i,Xj) is a function 
of the starting point Xi-i and ending point Xi of the i-th flight, and () denotes the average 
over all possible processes starting at point xq. Then, in analogy with the average number 
of flights, such a quantity must satisfy a recursion relation 

Q{x ) = (q {x )) + f Q(x 1 )p(x - x 1 )dx 1 , (13) 
Jo 

where (qo(x)) = f^ 30 p(xi — x )q(x ,Xi)dxi. Note that if x\ is outside the interval [0, L], 
the particle is absorbed by one of the boundaries and the value q(xo,x\) should be defined 
according to its physical meaning for the absorbed particle. Equation fll~3l ) is identical to 
Eq.fllTD with h(x) = -(q (x)). 

As an example of application of Eq. (|13"D , let us consider quantity Q(xq) to be the total 
flux through the right boundary. This flux is related to the transmission probability of 
photons through the clouds ||26|| . By definition, the flux through the right boundary is 



equal to the probability P r (xo) of the absorption of the flyer that starts at point x by the 
absorbing boundary x = L. In this case, quantity q must be defined as q(x , Xi) = 6{x\ —L). 
The very first flight is absorbed by the right boundary with probability p r {xo) = (qo(xo)) = 
/£° p{x\ — Xo)dxi, or after integration 

n(T) _j^o/(L-x )} a /2, 0<x o <L-£ o 

P^ x °)-\ 1/2, L-£ <x <L- li4j 



Therefore, P r (xo) satisfies Eq. flT3"| ) with (qo(x )) = p r (xo) and Eq. (JTTJ) with h(xo) = 
—p r (xo). In the next section, we will apply this method to define the total path length 
of the flyer. 

III. AVERAGE TOTAL PATH LENGTH TRAVELED BY THE FLYER 

The average total path length traveled by a Levy flyer before absorption is equivalent 
to the total time spent by a Levy walker before absorption [25||. The evolution of the 



probability density of Levy walks was studied |25j in terms of time. This approach leads to 
integral equations involving integration over time and space. Here we restrict our study to 
the problem of the average total path length before absorption of discrete Levy flights. This 
particular problem can be solved in much simpler terms. 

In the absence of the absorbing boundaries, the average flight length with probability 
density p(£) of Eq. (|j) is given by (\£\) = f_™ \£\p(£)d£, which is independent of the starting 
point and diverges for a < 1. In the presence of the absorbing boundaries, the flight starting 



6 



from a point y cannot exceed the distances y and L — y from this point to the boundaries. 
One can show that the average length of a flight that starts from a point y of an interval 
[£ ,L- £ ] is given by 

s(y) = (\t(y)\) = ^ 



dx 



y-to 

\y-xY 



+ 



dx 

y+to (x-y) a 



(15) 



dx 



00 dx 



,(y-x) a+1 v ' Jl (x-y) a+1 
which converges for a > 0. If a ^ 1, Eq. ([To] ) reduces to 



2(1 - a 

In case a = 1, we have from Eq. (IT 



a-l 



a-1 



,3/ J + U-?/, 



2a 



»(y) = J 



If < y < £ , Eq. ([TBI) must be replaced by 



/ x y . *o 

= o + 



ln(^ + ln^] + 2 



a-l 



or by 



2 2(1 -a) 



s(j/) = — — + 



L-y, 



a 



2 2fl -a) 



0-1 



y 



a 



(16) 



(17) 



(18) 



(19) 



if L — £ < y < L. Analogous changes must be made in Eq. ([Ti 

Thus, according to Eq. ([13]), the total average path length for the process that starts at 
point xq is 



or 



(S(x Q )) = s(x ) + / (S(x 1 ))p(x - xi)dx u 
Jo 



(S(x )) = -Uc a -l)- l s] (x ) 



(20) 



(21) 



This equation is identical to the Eq. (|H|) in which function h(x) = —s(x). We will solve 
Eqs. ( p!2|) and (|20|) numerically in Section V. 
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IV. THE CONTINUOUS LIMIT 



Appendix A shows that for a < 2, operator £ a {£o) — I tends to zero when £ —>■ for 
any function f(x) that has second derivative f"(x) for < x < L and has finite limits /(0) 
and f(L). It can also be shown that for such functions and a < 2 there exists an operator 

rag 

V a = lim V*[£ a (£ )-X]. (22) 

The result of this operator acting on any such function is defined as 
m f]( x vp /- L sgn(x - y)f'(x)dx f(Q) f(L) 

where /'(sc) is the first derivative of the function f(x). This operator is a self-conjugate 
operator similar to the double differentiation operator d 2 /dx 2 . It can be expressed |3,24 



as 



the linear combination of right and left Riemann-Liouville fractional derivatives of the order 
a. The difference of the two operators 

d a (£ )=£- a [C a (£ )-l]-V a (24) 

decays as £ 2 ] ~ a when £q — > 0. Appendix A also shows that the leading term of the operator 
d a is proportional to the operator of the second derivative 

In analogy with the diffusion equation with continuous time, we can define a superdiffu- 
sion equation [^|, |T2| , pL3|,|T6|l based on Levy flights. Instead of the discrete process defined by 



Eq. (pH, one can write 



dP{x,t) „ 



dt t 



-?-V a P(x,t), (26) 



where to is the duration of each flight, and ifi/to is the fractional analog of the diffusion 
coefficient. Note that t$ plays a role similar to the mean free path, and t plays the role of 
the mean collision interval. 

The operator T> a has an orthogonal normalized set of eigenf unctions fk{x), such that 
T^afk( x ) — ^^/^(x), and /fc(0) = fk(L) = |29|^0| . Similarly to the solution of usual 
diffusion equation, the solution of Eq. (E?j) can be expressed via separation of variables as a 
series of eigenfunctions 



00 , r L 

P(x,t) = Y,e X ^ t/to fk(x) f k (y)P(y,0)dy. (27) 
tZ Jo 



s W/* 0/fc(x) j L 

k=l 

In Ref. |]16|, where the method of separation of variables for the super diffusion equation 



on a finite interval has been first proposed, it has been assumed that the eigenvalues 
asymptotically behave at large k as ~ —(k/L) a < and that the eigenfunctions fk(x) 



S 



can be well approximated by the eigenfunctions y2/Lsm(xirk/L) of the Laplace operator 
with absorbing boundary conditions. Numerical studies |4l[] confirm these assumptions but 
show that eigenfunctions /& and sines have different behavior near absorbing boundaries, 
namely, fk(x) ~ x a ^ 2 as x — > 0. 

Having defined the properties of the operator T> a , we can derive the closed form expression 
for the average time spent by the continuous Levy flight process on the interval. Formal 
substitution of Eq. (|22|) into Eq. (O) yields 



(t)=to(n) = ^\v- 1 h] (* ) 



l 



where function g(x) satisfies the equation 



Dagix) = h(x) = -1. 



(28) 



(29) 



Note that g(x) has to satisfy boundary conditions g(0) = 0, g(L) = 0. Otherwise, according 
to Eq. ([23]), the right hand side of Eq. (^) would contain singularities. In the general case, 
the equation 



V p f L S ^(V - x)f'(y)dy _ 

Jo 2\y-x\ a V ' 



with absorbing boundary conditions 



/(0) = f(L) = 



(30) 



(31) 



belongs to a known class of generalized Abel integral equations with Riesz fractional kernel 
24,29-31 1 . It can be shown [29,31] that such an equation with boundary conditions ([31]) has 



a unique solution which can be obtained via spectral relationships for Jacobi polynomials 
24] , |3~ll| or by the Sonin inversion formula |]3T] , |32l| (see Appendix B). Similar inversion formulae 
are given in ref. [29]. In the case h(x) = —1, the solution can be expressed in elementary 
functions: 



g(x) = 2sin ^ a/2) [(L-x)x}^. 



(32) 



One can verify this solution by performing contour integration around the cut [0, L] on the 
complex plane and computing the residue of the integrand at infinity. It should be pointed 
out that g(x) can be expanded in a series of eigenfunctions fk{x): 



g{xo) 



oo , L 

Y.^fkixo) / fk(x)dx. 

k=l J ° 



This expansion is similar in spirit to an approximation, obtained in [03], where the exact 
eigenfunctions were approximated by sines. Although approximation [|16j correctly predicts 
the power law dependence g(xo) ~ L a for the points Xq in the center of the interval, it 
differs from Eq. ( |32|) in the proportionality coefficient and in the behavior near absorbing 
boundaries J|l] . 
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Note that for a > 2, Eq. (|26|) should be replaced (42| by the standard diffusion equation 
with diffusion coefficient D = £ 2 aj\t 2[a — 2)]. In this case, the average time spent by the 
flyer before absorption is given by the classical equation (t) = x(L — x)/{2D). Note that 
the diffusion coefficient D resembles the proportionality coefficient in Eq. (ffif). 

One can argue that the expression (|32| ) may yield the average number of flights taken by 
the discrete Levy flight process in the limit of £q — > 0. Indeed, according to Eq. (f2~4p 



where operator d a (£ ) 
obtain 



and thus 



£ a (£ )-I = £%[V a + d a (£o)), (33) 
0, as £ — > 0. Formally expanding (C Q — X) -1 in powers of d a , we 



; -(v- 1 -v- 1 d a v- 1 + ..: / 



(n)=£^{g(x )-[V- 1 d a g}(x ) + ...} 



(34) 



(35) 



Note that expansion fl3"5D is formal and may not converge. We will test this assumption 
numerically in Section V. In order to distinguish the average number of flights for the 
discrete process, (n), from the continuous limit approximation, we will denote the latter by 
n a (x): 



n a (x 



sin(7ra/2) 
vra/2 



[L - Xq)x 

f 2 



a/2 _ sm{na/2)M a {z - z 2 ) a l 2 
~ vra/2 : 



(36) 



where z = Xq/L, M = L/£ . 

Analogously, we will denote the continuous limit approximation for the average total 
path length by S a . In the continuous limit, Eq. ( pl~l) should be replaced by Eq. ( j30t ) with 
f(y) = S a (y), h(x) = £o a s(x), and absorbing boundary conditions (|3l|). In this case (See 



Appendix B), Sonin formula leads to an expression containing hypergeometric functions: 



S a (x ) 



L{2-a) 
2(1 -a) 



1-4 



ip a {z) + xpJl - z) 



'a a x 
• 2 ' 2 / 



2LM a ~ 1 sin ( ^ 



+ 



7r (a 



a(a + 2)B(^ 

where B(a, b) = T(a)T(b)/T(a + b) is Euler B-function, 
and F is the hypergeometric function |39| 

F(n + a)T(n + b)x r ' 



(37) 



a a a 
F(2- -,-,-+ 2, z)z 



2+1 
2 Ti 



F(a, b, c, x) = 



T(c) 



E 



r(o)r(6)^s r(n + i)r(n + c ) • 

In case a = 1, corresponding to the Cauchy distribution, the hypergeometric function ip a (z) 



can be expressed in elementary functions ipi(z) = | | + sin x (2z — 1) — 2V i — z 2 and 

Eq. (|37| ) yields S\{xq) = L-\fz — z 2 \nM + 0(1), where terms 0(1) do not depend on M 
and can be found from Eq. ([37]) by L'Hopital's rule. 
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Note that the average total path length traveled before absorption by the Levy flyer can 
be expressed in terms of survival probability 0(£) of the Levy walker: (S) = / °° Q(t)dt. 
According to J25[] , this survival probability exhibit for t — > oo asymptotic exponential decay 
0(t) ~ exp(— |Ai|t), where Ai is the main eigenvalue of the correspondent problem. Substi- 
tuting Q(t) by its asymptotic, we can estimate (S) ~ l/|Ai|. Thus (S) and l/|Ai| must have 
identical asymptotic behavior for large L. Indeed, Eq. (|37|) yields S a (x) ~ L for < a < 1 
and S a (x) ~ L° for 1 < a < 2, in complete agreement with asymptotic approximations of 



Finally, we will find the probability P r (xo) of the absorption by the right boundary in 
the continuous limit. According to Eqs. ( jl3]) and ([14]), P r (x ) should, in continuous case, 
satisfy Eq. (p9|) with h(xo) = ~£^ a p r (xo), i.e. 



V a T\ 



(38) 



For xq = 0, the flyer is immediately absorbed by the left boundary, so P r (0) = 0. For x = L, 
the flyer is immediately absorbed by the right boundary, so P r (L) = 1. Thus the second 
term in the expression (|23"D for I , a 'Pv(§/) is equal to zero and the third term cancels out with 
the right hand side. Hence P r (x ) satisfies homogeneous Eq. fl3U| ) and boundary conditions 
P r (0) = 0, P r (L) = 1. This solution can be expressed in terms of the homogeneous solution 
(fo(x) obtained in Appendix B: 



P r (x ) 



Io° <Po(y)dy 
Io ( Po(y)dy 



Xo\ 2 

LJ 



'a -| 
* 2 ' L 



a a 
2 ' 2 



+ 1 



XQ ' 

' L < 



a 

2 ' 2 > 



(39) 



Note that the probability of the absorption by the left boundary, Pi(xq) = P r (L — xq) = 
1 — P r (xo). For xq — > 0, the asymptotic behavior of P r (xo) is given by P r (xo) ~ {xq/ ' L) a l 2 
which is in complete agreement with the result of ref. |26| for the transmission probability 
of the photons through the clouds of depth L. 



V. NUMERICAL SOLUTION 



The goal of this section is to treat Eq. (12) and Eq. (20) numerically and to compare 
the results with the continuous limit solutions Eq. fl36|) and Eq. (|3TD. To perform numerical 
integration of Eqs.(|l2|) and (p0|) , we replace the integration by summation and the kernel 
p(x — y) by the matrix Ay, < i < M, such that Ai 



A 



Kl 



1 



j\ a (\i-j\ + iy 



and 
i ^ J- 



(40) 



Accordingly, the average flight length s(x) performed from the point k = x/£q is replaced 
by M — 1-dimensional vector s with elements 



Sk 



2(1 -a) 



k 



a-l 



+ 



(M -ky- 1 



- 2a 



(41) 



The average number of flights for the process that starts from point k = x/£o is 
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n) k = E A ^ ■ c = [(/ - A)- l c] u , (42) 



where c is the vector with all components equal to 1, and e& is a unit basis vector with 
components eu = Ski, where 5^ is the Kronecker delta. Analogously, average total length is 
equal to 




(I-A)^s\ k . (43) 

The symmetric matrix R = (A — J) -1 is the analog of the self-conjugate operator (C — X)~ l . 
Using iterative techniques for matrix inversion, we obtain the numerical solutions for (n) 
and (S). 

In Figure 1, we compare the numerical solution ( f42|) for (n(xo)) and the continuous 
limit approximation n a (xo) given by Eq. ( |36|) for the case xq = L/2. In order to test the 
asymptotic convergence, we have to divide both functions by (M/2) a . It can be seen that 
for Xo/£q ^> 1, Eq. ([36]) provides good approximation for the average number of flights of the 
discrete process defined by Eq. ([TTJ) . Studying the difference between the numerical values of 
(n)(M/2)~ a and the continuous approximation 2 sin(7ra/2)/(7ra) for xq = L/2, we confirm 
that for 2 > a > 1 this difference decays as M a ~ 2 with M — > oo. This is in agreement 
with Eq. (p5|). However, for a < 1 the difference between the numerical solution and the 
continuous approximation converges as M _1 > M a ~ 2 . The term M _1 is proportional to the 
error of replacement of integration in Eq. (|1^) by summation. 

In Fig. 2, we compare the numerical solution ( f4~3"D for (S(xq)) and the continuous limit 
approximation S a (xo) given by Eq. (|37| ) for the case xq = L/2. In order to test the 
asymptotic convergence in this case, we have to divide both functions by (LM Q_1 —L)/{a — 
1). It can be seen that for x /£q 3> 1, Eq. Q3TD provides good approximation for the average 



total traveled length in the discrete process defined by Eq. (|2T|) 

Now we will examine the quality of the continuous limit approximation in the vicinity of 
the absorbing boundary. For simplicity we will study only the behavior of the average number 
of flights. The approximation for the total path length has similar problems. Figure 3 shows 
that for Xq/£q = 1, the correction terms in Eq. ( pop cannot be neglected. As shown in 
Appendix A, the operator d a — > for any fixed x if £ — > 0, but does not vanish if x and £ 
both approach zero, so that their ratio x/£q — > r > 0. Accordingly, the value of (n) behaves 
for xo/£o = r, M — > oo as x(a)(rM) a//2 , where x{°) > 2sin(7ra/2)/(7ra) is some unknown 
function that can be estimated numerically (see Fig. 3). It is likely that x( a ) remains 
positive as a — ► 2. The analytical determination of the function x( a ) remains an unsolved 
problem. Nevertheless, continuous approximation correctly predicts the leading factor M Q//2 
for the average number of flights started in the vicinity of the absorbing boundary. 

In summary, comparison of the numerical solutions for the mean first passage time of 
Levy flights and Levy walks and the exact solutions of these problems in the continuous 
case suggests that fractional differential equation for super diffusion Eq. ( p6|) with absorbing 
boundary conditions provides good approximation for discrete Levy flights on a finite interval 
with absorbing boundaries. However this approximation breaks down when a — ► 2 and in 
the vicinity of the absorbing boundaries. 
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VI. ANALYSIS OF THE TOTAL PATH LENGTH: IMPLICATIONS FOR 

BIOLOGICAL FORAGING 



Recently, biological foraging has been modeled by Levy flights [|7|,|33|-p^] . The case of 
non-destructive foraging (defined in |36| as case in which "target sites" can be revisited not 
just once but many times) corresponds to x = £ , i.e. the forager starts its next search 
from the previously-visited food site, located at the origin. The prey may reappear at this 
site. Accordingly, coming back to the origin may be profitable in terms of foraging efficiency, 
which is defined [^6| as the inverse average total path length before finding next food site. 
With help of Monte-Carlo simulations, it has been shown [[SB] that, in case of non-destructive 
foraging, the foraging efficiency has maximum at a = 1. 

We confirm this result, using numerical solution ([43"1). Figure 4 shows a semi-logarithmic 
plot of (S) versus a for xq = £q and various values of M. On can see that the minima, 
a m i n (M), shift towards a = 1 as M — ► oo. Heuristic approximations |36|,|37j suggest that 
foraging efficiency has a maximum at 1 — e, where e ~ (lnM)~ 2 . Consequently, the average 
total path length should have minimum at the same point. Figure 5 confirms this prediction 
for the numerical solution. It shows the graph of ct min (M) versus [ln(M)]~ 2 , which is almost 
a straight line with an intercept a w^ oo) ~ 1. 

In the following, we will prove this result using the continuous limit approximation flT7|). 
Accordingly, we will find the behavior of S a (xo) for the case when the starting point xq is 
selected in the vicinity of the absorbing boundary and show that in this case S a (xo) has a 
minimum at a — > 1. In oder to do this, we will present solution ( |37| ) in a more convenient 
form, which allows to separate leading singularities at z — ■> 0. After some transformations 
involving hypergeo metric functions [39], we can rewrite Eq. (37) as follows: 



S a (xo) 



2L(z 



a 



2\ — 

Z ) 2 



M a-l gin / t„ 



0) {l-z)h{^z) 



7T 



oB(f,f) 



(44) 



L(2 
2(cT 



a) 



ft vi 4/ 3 (a,2;)z2 
f 2 {a,z)(l- z)* + - c 



where 



OL 

f 1 {a,z) = F(a,2,- + l,z), 



and 



f 2 (a,z) = F(2-^,l-^,z), 



^/ a a a 
f a ( a ,z) = F(2 --,-,- + 2, z). 



The first term in solution (^3) decreases as z a ^ 2 , when z — * 0, while the second part decreases 
as z. If we take the starting point xq = r£o, where r is constant, then z = r/M, where M is 
large number, and we can separate the leading (with respect to M) terms in solution 
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S a (r£ ) = r Q / 2 Ua)M%- 1 - ((a)M~%] + L ■ OfM" 1 ), (45) 

a — 1 L - 1 

where 

. , 2sin(-7ra/2) 2 . . 

V^) = C(«) = -575 ay ( 46 ) 

The above approximation accurately follows the solutions of the discrete problem for a < 1, 
when the term £(a)M _Q//2 dominates, but strongly deviate from that for a > 1, when the 
term r](a)M^~ l dominates (see Fig. 4). The reason for these deviations is the truncation 
of the non-leading terms in the Eq. 



In contrast with the discrete solution and Monte-Carlo simulations of [36], expression 
Eq. (HSp has two minima: one at a = 1 + e(M) and another at a = 2. We will show that 
e(M) — > as M — > oo. Let us expand r\ and £ in powers of e: 



??(! + e) = r/o + + C(l + e) = Co + &e + 



(47) 



Note that 7/ = Co = 2/-7T, and hence the expression ( fR-j ) does not have a singularity at a = 1. 
The location of the minimum can be found by differentiation of the expansion for Si +e (r£o) 
with respect to s and equating the leading terms of the order of In M: 

Indeed, Eq. ([48]) shows that e(M) — > as M — > oo and r stays constant. 

This analysis holds for any continuous functions 77 and (, so long as 7/(1) = ((1) and, 
therefore, is likely to be valid in the discrete case, in which, functions 77(a) = rjd(a) and 
C(a) = Cd(&) do not satisfy Eq. (|46|) . Note that r] d (a) can be expressed in terms of function 
x(a) shown in Fig 3., namely r] d (a) = ax{cn). Analysis of Fig 3. shows that x(2) > 0. 
Consequently, the minimum at a = 2 does not exist in the discrete case. 



VII. SUMMARY 

We have studied Levy Flights in a finite interval with absorbing boundaries. In Section 
II, we have derived expressions Eqs. ([TTJ) and ( VI ) for the average total number of flights 



(mean first passage time). We also obtain a general recursion relation Eq. ( |T3"D for the 
average of the sum of arbitrary contributions from each flight in the form of the Fredholm 
integral equation of the second kind. We applied this method to derive the probability of 
absorption by one of the boundaries. In section III, we have derived expressions Eqs. ( |20|) 
and (El]) for the average total path length of the Levy flights which is equivalent to the mean 
first passage time of the Levy walks. 

In Section IV, we have shown (See Appendix A) how the discrete Levy flights are related 
to the fractional differential equation Eq. (|26D of the superdiffusion with Riesz operator 
Eq. fl2"5p. For the continuous process described by Eq. (EH), we derived exact analytical 
expressions Eqs. (|36D, fl3~7]) , and d39|) for the mean first passage time, the average total path 
length, and the probability of absorption by one of the boundaries, respectively. All these 
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quantities are the solutions (See Appendix B) of the fractional differential equation Eq. fl3"0| ) 
with Riesz kernel and with different right hand sides. In Section V, we have compared these 
analytical solutions with numerical solutions obtained for the discrete Levy flights (See Figs. 
1-3). We have shown that fractional differential formalism provides good approximation for 
the discrete Levy flights in the interval with absorbing boundaries except the case of a —>■ 2 
and the case when the starting point is in the vicinity of an absorbing boundary. In the 
latter case the fractional differential formalism yields correct scaling behavior with respect 
to the interval size and distance to the boundary, but gives an incorrect proportionality 
coefficient (See Fig. 3). 

In Section VI, we have investigated the behavior of the average path length as a function 
of the starting point and as a function of a. We have derived asymptotic expression Eq. 
(|45|) for this quantity in the case when the starting point is located close to the absorbing 
boundary. We have shown that the expression for the average path length has a minimum at 
a ~ 1 if the process starts in the vicinity of the absorbing boundaries (See Figs. 4,5). This 
result, as well as Eqs. (|39|) and (f44|), can be applied to the problem of light transmission 
through cloudy atmosphere p6|j27|] . 

Similar fractional integral operators |3|,|T3"|1 - namely the Riesz operator V a — can be 
used to treat the problem of the Levy flyer in the dimensions higher than one with randomly 
distributed absorbing traps. Let L be a characteristic distance between neighboring traps. 
Then we still expect that the average number of flights before absorption scales as L a if the 
process starts far away from the absorbing boundary and as L a l 2 if the process starts in the 
vicinity of the absorbing boundary. This result is sufficient to prove that the minimum of 
the average total path length traveled by the flyer before absorption is achieved at a — > 1, 
if the flyer starts in the vicinity of the absorbing point. 

Finally, we comment on the relevance of our findings to biological Levy flight foraging. 
Our results essentially confirm that Levy flights with a = 1 (or fi = 2 in notation of 
should theoretically provide the optimal strategy of foraging in case of sparsely and 



randomly located food sites, if any food site can be revisited many times ||36|| . The presence of 
the second minimum near a = 2 predicted by continuous limit approximation may indicate 
another possible strategy for foraging, i.e to perform Brownian walks in the region of possible 
appearance of prey. Break down of the continuous limit approximation in the vicinity of the 
absorbing boundary indicates that the results should depend on the particular details of the 
model. 
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APPENDIX A: EXISTENCE OF THE CONTINUOUS LIMIT OPERATOR 



We will show that operator D a is defined on any function f(x) that has finite limits at 
both ends of the interval /(0) and f(L) and finite second derivative f"(x) at any inner point 
x of the interval [0, L\. According to Eqs. (|]) and (|22|) , 



Vafiv) = Km /o 



f(x)dx 



y- £ o f(x)dx r L 

(y - x) a+l + iy+fo (x - y) Q+1 



- f(v) ■ 



(Al) 



Making partial integration of both integrals in Eq. (|A1|) we get 

to— »0 



-/(v) 



+ - 



f-^o f'(x)dx r L f(x)dx 
(y-x) a J y +e (x-y) a 



(A2) 



/(0) /(L) 



2(L - y) c 

For a < 2, the first term in Eq. (|A2[) goes to zero as 
The second term converges to the integral 



I = V.P. 



sgn(a; — y)f'(x)dx 



2\y 



x 



(A3) 



which exists for a < 2 if /'(x) has a derivative at a: = y. 

Subtracting Eq. ( |A3|) from the second term in Eq. (|A"2j), replacing i n ^he integrand 
by its Taylor expansion f'(x) = f'(y) + f"{y){x — y) + o(x — y), and combining it with the 
first term of Eq. ( |A2| ) , we reproduce Eq. for the correction operator d a (£o). This shows 
that the operator D a is well defined for the class of functions with existing second derivative. 



APPENDIX B: SONIN INVERSION FORMULA 

Equation (|30|) belongs to a class of generalized Abel equations. In his classical works, N. 
Ya. Sonin [|EJ suggested a general method for solving such equations. In particular |31 |, an 
equation: 



L [aisgn (x-y) + a 2 ] 
2 | x - y \ a 



-ip{y)dy = h(x) 



< x < L 



(Bl) 



has a solution: 
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A r L fl-ot A r t n -/ 3 

^ = B --V z I ^7r-i t jt I v^ h ^ (B2) 

where 

B a = -2 sin (vr/5) T (a) r^r -1 (1 - (3 + a) (cii + a 2 y v tt -1 , (B3) 
and parameter (3 is determined by relations 

sin (7r(/5 — a)) = c u sin {irf3) , (B4) 



(12 — (l\ 

a\ + a 2 



(B5) 



Similar inversion formulae can be found in ref. |29| . 

In case of Eq. (p0|), v — l,a2 = 0, a\ — —1. Hence, according to Eq. (p5| ) c v = — 1. 
Equation (|B4j ) has infinite number of solutions (3 = a/2 + k, where k is an integer. For 



< a < 2, only two solutions with k — 0, k — 1 lead to the converging integrals in Eq. (|B2] 

r/ r L , d r* 



and 



BaZa ' 2 ~dzl z dt^it-zy^jtly^it-yT^hWy ( B6 ) 



WJ = ~B a z^±- £ dt t l -(t - zY /2 - l jjl V a/2 -\t - y^ 2 h(y)dy, (B7) 



where 



4sin(^) , , 



vraB(f 



2 ' 2 • 



Since Eq. contains /'(x) = v 9 ^); one can always satisfy the first boundary condition 
Eq. (j3l|) by defining f(x) = J X ip(z)dz. Adding solution (|B6|) for h(x) = —1 and solution 
( |B"7j ) for /i(a;) = 1, one can see that the homogeneous equation with h(x) = has a nontrivial 
solution tpo = (Lx — x 2 )? -1 . Hence the second boundary condition Eq. ([H]) can be satisfied 



if we select <p = tpi — Cip , with constant C = L l a J L (p(x)dx/B(^ } |). In case h(x 



straightforward calculations lead to Eq. fl32|). 

Now we will obtain the analytical solution for the average total path length before ab- 
sorption in the continuous process, S a (x). In this case, the right hand side of Eq. (|30|). 
h(x) = s(x), is given by Eq. (|16|), and we can use its symmetry h(x) = h{L — x). Thus, in 
order to satisfy the second boundary condition Eq. ([H|), we should have <p(x) = — f(L — x). 
To construct such a solution, we first find the solution ipi(x) for the first term in Eq. ([16]) 
h(x) = Obviously, the function —ipi(L — x) provides the solution for the second term 

h(x) = (L — x) 1_Q . The solution for the third constant term is given by Eq. ( p2[) with a 
proper coefficient. Summing up all three partial solutions and using various properties of 
hypergeometric functions p9| , one can find the total solution presented in Eq. ([37]). 
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FIG. 1. The behavior of the scaled average number of flights (n)(M/2) a versus a for in- 
creasing values of M = 50, 100,200,400,800, 1600,3200 in the case xq = L/2 in comparison with 



the continuous limit prediction of Eq. (36), 2sin(7ra/2)/(7ra), shown as a bold line. We see good 
convergence to the predicted function except for the values of a ~ 2. We extrapolate the values 
(n)(M/2)~ a for M — > oo (circles) using their polynomial fits with respect to M _1 for a < 1, or 
with respect to M 2 ~ a for 1 < a < 2, or with respect to (InM) -1 for a = 2. 
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FIG. 2. The behavior of the scaled average total path length (S)(a — 1) / (LM a ~~ 1 — L) versus 
a for increasing values of M = 50, 100, 200, 400, 800, 1600 in the case xq = L/2 in comparison with 
the continuous limit prediction of Eq. ( |37| ) in the limit M — ► oo shown as a bold line. In this limit, 
continuous approximation follows first term of Eq. (]37| ) for a < 1 and the second term sin(7ra/2)/7T 
for a > 1. We see good convergence to the predicted function except for the values of a « 2. 
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FIG. 3. The behavior of the scaled average number of flights (n)(M)~ a / 2 versus a for in- 
creasing values of M = 50, 100, 200, 400, 800, 1600, 3200 in the case xq = £o in comparison with 
the continuous limit prediction of Eq. ( |36| ) 2 sin(7ra/2)/(7ra), shown as a bold line. Although the 
values are close to the continuous limit predictions, they converge to a different function x( a ) as 
M — * do. To obtain x( a )> we extrapolate the values {n)(M/2)~ a for M — * oo using the same 
procedure as in Fig. 1. We assume that the error bars are equal the discrepancies between the 
extrapolation and the continuous limit in Fig. 1. 
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FIG. 4. Semi- logarithmic plot of numerical solutions of (S) versus a for the case xq = £q = 1 
and various values of M = 50, 100, 200, 400, 800. Circles indicate the positions of the minima 
Oimin(M) which shift towards the vertical line a = 1, as M increases. In addition, we show the 
analytical continuous limit approximation S a (£o) given by Eq. (46) for M = 800. 
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